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High order finite difference methods with subcell 
resolution for 2D detonation waves f 

By W. Wang}, C.-W. Shut, H. C. Yee|| AND B. Sjogreenff 


1. Motivation and objective 

In simulating hyperbolic conservation laws in conjunction with an inhomogeneous stiff 
source term, if the solution is discontinuous, spurious numerical results may be produced 
due to different time scales of the transport part and the source term. This numerical 
issue often arises in combustion and high speed chemical reacting flows. 

The reactive Euler equations in two dimensions have the form 

U t + F(U) x + G(U) y =S(U), (1.1) 

where U, F(U), G(U) and S(U) are vectors. If the time scale of the ordinary differential 
equation (ODE) Ut = S(U) for the source term is orders of magnitude smaller than 
the time scale of the homogeneous conservation law Ut + F(U) X + G(U) y = 0 then 
the problem is said to be stiff. In high speed chemical reacting flows, the source term 
represents the chemical reactions which may be much faster than the gas flow. This leads 
to problems of numerical stiffness. Insufficient spatial/temporal resolution may cause 
an incorrect propagation speed of discontinuities and nonphysical states for standard 
dissipative numerical methods. 

This numerical phenomenon was first observed by Colella et al. (1986). Then LeVeque 
& Yee (1990) showed that a similar spurious propagation phenomenon can be observed 
even with scalar equations. Colella et al. (1986) and Majda & Roytburcl (1990) have 
successfully applied the random choice method of Chorin (1976, 1977) for the solution 
of under-resolved detonation waves. However, it is difficult to eliminate completely the 
numerical viscosity in a shock-capturing scheme. Fractional step methods are commonly 
used for allowing an under-resolved mesh size with a shock-capturing method. Chang 
(1989, 1991) applied the subcell resolution method of Harten (1989) to the finite volume 
ENO method in the convection step, which is able to produce a zero viscosity shock profile 
in nonreacting flow. The time evolution is advanced along the characteristic line. Correct 
discontinuity speed was obtained in the one-dimensional scalar case. However, it is diffi- 
cult to extend this approach to multi-dimensions and to system of equations because of 
the reliance on the exact time evolution via characteristics. Engquist & Sjogreen (1991) 
proposed a simple temperature extrapolation method based on finite difference ENO 
schemes with implicit Runge-Kutta time discretization, which uses a first- /second-order 
extrapolation of the temperatures from outside the shock profile. The method is easy to 
extend to multi-dimensions. Their method is not a fractional step method. It does not 

f This is a condensed version of “High order finite difference methods with subcell resolution 
for advection equations with stiff source terms” in J. Comput. Phys., to appear. 
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seem to work well when the spatial scales are under-resolved. Other first-/second-order 
methods that are based on the fractional step method have been proposed by Bao & Jin 
(2000, 2001) and Tosatto & Vigevano (2008). 

In our previous work Wang et al. (2010), we developed a high-order finite difference 
method which can capture the correct detonation speed in an under-resolved mesh and 
will maintain high-order accuracy in the smooth part of the flow. Numerical examples 
were presented for one dimensional scalar problems and one dimensional detonation 
waves. Our objective in this study is to extend this method to two dimensional reac- 
tive Euler equation. The first step of the proposed fractional step method is the con- 
vection step which solves the homogeneous hyperbolic conservation law in which any 
lrigh-resolution shock-capturing method can be used. The aim in this step is to produce 
a sharp wave front, but some numerical dissipation is allowed. The second step is the 
reaction step where an ODE solver is applied with modified transition points. Here, by 
transition points, we refer to the smeared numerical solution in the shock region, which is 
due to the dissipativity of a shock-capturing scheme. Because the transition points in the 
convection step will result in large erroneous values of the source term if the source term 
is stiff, we first identify these points and then extrapolate them by a reconstructed poly- 
nomial using the idea of Harten’s subcell resolution method. Unlike Chang’s approach, 
we apply Harten’s subcell resolution in the reaction step. Thus our approach is flexible 
in allowing any shock-capturing scheme as the convection operator. In the reaction step, 
since the extrapolation is based on the high order reconstruction, high order accuracy can 
be achieved in space. The only drawback in our current approach is that the temporal 
accuracy will only be, at most, second-order due to the time splitting, which is common 
for most of the previous methods for stiff sources. We also remark that, in order to resolve 
the sharp reaction zone, sufficiently many grid points in this zone are still needed. The 
proposed method can capture the correct location and jump size of the reaction front, 
but it does not resolve the narrow reaction zone as typically there is one or few point in 
that zone. 


2. Two-dimensional reactive Euler equations 

The considered two-dimensional problem is modeling the reaction with two chemical 
states: burnt gas and unburnt gas. The unburnt gas is converted to burnt gas via a single 
irreversible reaction. Without heat conduction and viscosity, the system can be written 
as 


Pt + ( pu) x + ( pv)y = 0 

(2.1) 

(, pu)t + ( pu 2 +p) x + ( puv)y = 0 

(2.2) 

(, pv)t + (puv) x + (pv 2 +p)y = 0 

(2.3) 

Et + ( u(E + p)) x + ( v(E + p)) y = 0 

(2.4) 

(pz) t + ( puz) x + ( pvz) y = —K(T)pz, 

(2.5) 


where p(x,y,t ) is the mixture density, u(x,y,t) and v(x,y,t) are the mixture x - and y- 
velocities, E{ x, y , t) is the mixture total energy per unit volume, p(x, y , t) is the pressure, 
z(x,y,t) is the mass fraction of the unburnt gas, K(T ) is the chemical reaction rate and 
T{x,y 1 t) is the temperature. The pressure is given by 

P=( 7 - 1 )(E - \p(.u 2 + v 2 ) - q 0 pz ), 


( 2 . 6 ) 
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where the temperature T = 2 and qo is the chemical heat released in the reaction. 

The reaction rate K{T) is modeled by an Arrhenius law 

K(T) = K 0 e X p(^^), (2.7) 


where Kq is the reaction rate constant and T) gn is the ignition temperature. The reaction 
rate may be also modeled in the Heaviside form 


K(T) 


1 /e T > Ti gn 
0 T < Ti gn ’ 


(2.8) 


where e is the reaction time and 1/e is roughly equal to Kq. For simplicity, we only 
consider the Heaviside source term (2.8). 


3. Numerical method for two-dimensional reactive Euler equations 

The general fractional step approach based on Strang-splitting (Strang 1968) for equa- 
tion 

U t +F(U) x + G(U) v =S(U), (3.1) 

is as follows. The numerical solution at time level t n+ i is approximated by 


U”+i = A (^p) R(At)A 


(3.2) 


The convection operator A is defined to approximate the solution of the homogeneous 
part of the problem on the time interval, i.e., 


U t + F{U) x + G(U) y = 0, t n <t<t n+1 . (3.3) 

The reaction operator R is defined to approximate the solution on a time step of the 
reaction problem: 

d 4 = S(U), t n <t<t n+1 . (3.4) 

at 

The convection operator is over a time step At and the reaction operator is over At/2. 
The two half-step reaction operations over adjacent time steps can be combined to save 
cost. 

Next, we introduce the proposed fractional step methods for the convection step and 
the reaction step separately. 


3.1. Convection operator 

In the convection step, we use fifth-order WENO (WEN05) with Lax-Friedrichs flux and 
third-order Runge Kutta for time discretization. 

3.2. Reaction operator 

If there is no smearing of discontinuities in the convection step, any ODE solver can be 
used as the reaction operator. However, all the standard shock-capturing schemes will 
produce a few transition points in the shock when solving the convection equation. These 
transition points are usually responsible for causing incorrect numerical results in the stiff 
case. Thus we cannot directly apply a standard ODE solver at these transition points. 
Here we use Harten’s subcell resolution technique in the reaction step. The general 
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idea is as follows. If a point is considered a transition point of the shock, information 
from its neighboring points which are deemed not transition points will be used instead. 

In the two-dimensional case, we apply the subcell resolution procedure dimension by 
dimension. 

(1) Use a “shock indicator” to identify cells in which discontinuities are believed to 
be situated. We consider the minmod-based shock indicator in Harten (1989); Shu & 
Osher (1989). Identify troubled cell / i; in both x- and y-directions by applying the shock 
indicator to the mass fraction z. 

Define the cell Ijj as troubled in the ^-direction if |s-| > and \sfj\ > |sf +lj - 

with at least one strict inequality where 

Sij = minmod{wj + ij - Zij,Uij - Zi-ijj- (3.5) 

Similarly we can define the cell I l3 as troubled in the y-direction if | s v - | > -_ 1 | and 

|«V-| > \ s y ij+1 \ with at least one strict inequality where 

= minmod{ 2 iiJ+ i - Zij, z i3 - Zij-ij- (3.6) 

If lij is only troubled in one direction, we apply the subcell resolution along this 
direction. If I t j is troubled in both directions, we choose the direction which has a larger 
jump. Namely, if \sfj\ > \ s\j | , subcell resolution is applied along the ^-direction, otherwise 
it is done along the y-direction. 

In the following steps (2)-(3), without loss of generality, we assume the subcell resolu- 
tion is applied in the x-direction. 

Assuming I l3 is troubled in the x-direction, we apply subcell resolution along the x- 
direction. 

(2) In a troubled cell identified above, we continue to identify its neighboring cells. 
For example, we can define Ii+ip as troubled if |s* +1 -| > l s f-i.jl and |sf +1 -| > | sf +2 -| 
and similarly define Ii-ij as troubled if |s*_! -| > |sf_ 2 -| and \sf_ 3 -| > |sf +1 -|. If the 
cell Ii- S ,j and the cell I t + r ,j ( s , r > 0) are the first good cells from the left and the right 
(i.e., Ii- s +ij and J i+r _ ij are still troubled cells), we compute the fifth order ENO in- 
terpolation polynomial Pi- S j(x) and Pi+ r ,j(x) for the cells and h+ r ,j , respectively. 
Because of the anti-diffusive corrector in the convection step, r and s will not be larger 
than 2 in general. 

Modify the point values Zij, T,j and pij in the troubled cell I t j by the ENO interpolation 
polynomials 

/ Aj — l^ij Pi— s,j (^i ) 'l '') ? Pij — Pi—s y j{p^iip)i if ^ ^ ^ 

\ ^ij — Pi+rj{Xi,z), Tij Pi-\-r.'i (x v ;, l ) : Pij Pi-\-r.'i (x,; , ()) . if 0 <C X v ; 

where the location 9 is determined by the conservation of energy E 

r@ r x i+ 1/2 

/ pi- s j(x\ E)dx + / pi +r j(x\ E)dx = Eij Ax. (3.8) 

J Xi— 1/2 ^ & 

Under certain conditions, it can be shown that there is a unique 9 satisfying Eq. (3.8), 
which can be solved using, for example, a Newton’s method. If there is no solution for 9 
or there are more than one solution, we choose Zij = Zi+ r j, T = Ti +r j and/O^ = Pi+ r ,j- 
However, in the system case we would like to have the shock travel ahead of the reaction 
zone, so we take the values of 2 , T and p ahead of the shock. 



5 


High order methods for 2D detonation waves 
For simplicity, in the considered stiff problem, the value of Z{ 3 can be taken as 


f 0, 6 > Xi 

\ 1 , 9 < Xi ' 


(3.9) 


(3) Use Uij instead of U t j in the ODE solver if the cell f,j is a troubled cell. 

For simplicity, explicit Euler is used as the ODE solver. 

(pz)Z +1 - ipz) Z + A tSff^pij, Zij). (3.10) 

Here we would like to remark that, implicit methods cannot be used in this step 
because the troubled values need to be modified explicitly. However, there is no small 
time step restriction in the explicit method used here, because once the stiff points have 
been modified, the modified source term S (Tij , pij , z^j) is no longer stiff. Therefore, a 
regular CFL number is allowed in the explicit method. 

In general, a regular CFL=0.1 can be used in the proposed scheme to produce a stable 
solution. But the solution is very coarse in the reaction zone because of the underresolved 
mesh in time. In order to obtain more accurate results in the reaction zone, we evolve 
one reaction step via N r sub steps, i.e. , 



in some numerical examples. 


3.3. A Numerical example of 2D detonation waves 
This example is taken from Bao & Jin (2000). The chemical reaction is modeled by the 
Heaviside form with the parameters 

7 = 1.4, q 0 = 0.5196 x 10 10 , ^ = 0.5825 x 10 10 , T ign = 0.1155 x 10 10 


in CGS units. 

Consider a two-dimensional channel of width 0.005, the upper and lower boundaries 
are solid walls. The computational domain is [0, 0.025] x [0, 0.005]. The initial conditions 
are 


{p,u, v,p, z) 


(Pl,ui, 0,p;,0), iix<f(y), 

(pr,u r ,0,p r ,l), ifx>£(y), 


(3.12) 


where 


t(y) = 


0.004 \y- 0.0025| <0.001, 

0.005 - | y- 0.0025| j y - 0.0025| < 0.001, 


(3.13) 


and ui = 8.162 x 10 4 , pi = 1.201 x 10“ 3 and Pl = 8.321 x 10 5 . 

Similar problems are also computed in Engquist & Sjogreen (1991). One important 
feature of this solution is the appearance of triple points, which travel in the transverse 
direction and reflect from the upper and lower walls. A discussion of the mechanisms 
driving this solution is given in Kailasanath et al. (1985). 

Figures 1-2 show density contours computed by WEN05/SR with 500 x 100 (Ax = 
Ay = 5 x 10“ 5 ), CFL=0.1 and N r = 2 at eighteen evolutionary times from t = 0 
to t = 1.7 x 10“'. We can see the movement of the triple points. The same case by 
WEN05/SR with a much coarser grid 200 x 40 (Ax = Ay = 1.25 x 10“ 4 ) with CFL=0.1 
and N r = 2 at three evolutionary times are shown in Fig. 3. We can see WEN05/SR 
with the very coarse 200 x 40 mesh can still capture the correct shock location, although 
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the shocks are smeared due to the lack of resolution. It is more apparent to compare the 
computed results with the reference solution in a ID cross section. The reference solutions 
are computed by standard WEN05 with 2000 x 400 grid points and CFL=0.3. The results 
by WEN05/SR and the splitting WEN05 are compared with the same mesh 200 x 40 
and CFL=0.005. Figures 4-7 show the ID cross section at y = 0.0025 at evolutionary 
times t = 2 x 10 -8 , t = 6 x 10 -8 , t = 1.4 x 10~ 7 and t = 1.7 x 10“ ' separately, where 
the left subplots are computed by WEN05/SR and the right subplots are by splitting 
WEN05. We can see WEN05/SR has excellent agreement with the reference solutions 
except it cannot capture the waves sharply due to the underresolved mesh. However 
the splitting WEN05 method produces spurious waves in front of the detonation shock 
starting at time t = 2 x 10 -8 (right subplot of Fig. 4) and after that the solutions move 
at a wrong speed (right subplots of Figs. 5-7). 


4. Future plans 

In this report, we demonstrated that the proposed high-order finite difference schemes 
with subcell resolution are able to capture the correct discontinuity speed in both one- 
dimensional scalar and system cases. Future work will extend this approach to higher 
dimensions with more general chemical reaction models. We will also consider multiple 
reaction models. 
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Figure 2. Computed density: WEN05/SR with 500x 100, CFL=0.1 and N r = 2 at nine different 
evolutionary times t = 9 x 10 s , t = lx 10 7 , t = 1.1 x 10 — 7 , t = 1.2 x 10 7 , t = 1.3 x 10”', 
t = 1.4 x 10” 7 , t = 1.5 x 10~ 7 , t = 1.6 x 10 -7 and t = 1.7 x 10 -7 . 
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Figure 3. Computed results: WEN05/SR with 200 x 40, CFL=0.1 and N r = 2 at 
t = 1.5 x 10~ 7 , t = 1.6 x 10~ 7 and t = 1.7 x 10~ 7 . 




Figure 4. ID cross-section at t = 2 x 10 -8 by different WENO schemes with 200 x 40. Solid line: 
reference solution; dashed line: numerical solution. Left: WEN05/SR with CFL=0.1, N r = 2; 
right: splitting WEN05 with CFL=0.05. 




Figure 5. ID cross-section at t = 6 x 10 -8 by different WENO schemes with 200 x 40. Solid line: 
reference solution; dashed line: numerical solution. Left: WEN05/SR with CFL=0.1, Ay = 2; 
right: splitting WEN05 with CFL=0.05. 
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Figure 6. ID cross-section at t = 1.4 x 10“' by different WENO schemes with 200 x 40. 
Solid line: reference solution; dashed line: numerical solution. Left: WEN05/SR with CFL=0.1, 
N r = 2; right: splitting WEN05 with CFL=0.05. 




Figure 7. ID cross-section at t — 1.7 x 10 'by different WENO schemes with 200 x 40. 

Solid line: reference solution; dashed line: numerical solution. Left: WEN05/SR with CFL=0.1, 

N r = 2; right: splitting WEN 05 with CFL=0.05. 
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